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Ferropericlase, (Mg,Fe)0, is the second-most abundant mineral of the Earth’s lower mantle. With 
increasing pressure, the Fe ions in the material begin to collapse from a magnetic to non-magnetic 
spin state. We present a finite-temperature first-principles phase diagram of this spin crossover, 
finding a broad pressure range with coexisting magnetic and non-magnetic ions due to favorable 
enthalpy of mixing of the two. Furthermore, we find the electrical conductivity of the mineral to 
reach semi-metallic values inside the Earth. 


Ferropericlase, (Mgi_a;Fea;)0, is an Fe-bearing transi¬ 
tion metal oxide that makes up some 20% of the total 
volume of the Earth’s lower mantle [H. Each Fe ion in 
this mineral assumes an octahedral coordination environ¬ 
ment, which leads to crystal field splitting, i.e., separa¬ 
tion of the Fe 3d shell of electrons into a higher-energy 
eg and a lower-energy t 2 g group. At low pressure, the 
ground-state electronic configuration of Fe^'*' is a high- 
spin state with four unpaired electrons giving a total 
spin of S' = 2. On compression, three effects come into 
play which ultimately cause a spin transition or mag¬ 
netic collapse to the low-spin state, S = 0. Firstly, the 
crystal field splitting grows due to increased overlap of 
the Fe and O valence orbitals, while secondly, the elec¬ 
tronic bands are broadened in energy due to increased 
confinement, making the high-spin state increasingly un¬ 
favorable ii- Thirdly, the low-spin state is favored 
by the smaller size of the low-spin Fe ion via the PAV 
contribution to the free energy Hi- 


Ever since the discovery of the spin transition in 
Feo .940 beyond 60 GPa at room temperature 0 and later 
in ferropericlase between pressures of 50 to 70 GPa 0, 
it has become a ppa rent that the phenomenon affects 
mechanical Is -111, compositional Q , and electronic 
properties [12l . Il3 |. and thus holds potentially signifi¬ 
cant implications for the physics and chemistry of the 
Earth. Interest in spin transitions is not, however, lim¬ 
ited to geoscience, with applications in, e.g., nanoclusters 
and thin films showing great technological potential 141. 
Experimental work utilizing x-ray emission spectroscopy 
(XES) Hi, optical spectroscopy [13: Mdssbauer spec¬ 
troscopy (MSB) 0, and equation-of-state (EOS) data 
gathered from high-pressure x-ray diffraction experi¬ 
ments [l^ has, to date, probed the spin transition in fer¬ 
ropericlase up to pressures of P = 140 GPa and tempera¬ 
tures of T = 2000 K. On the theoretical side, approaches 
based on analytical mean-field theory [13 and static first- 
principles calculations ll^ augmented by quasi-harmonic 
phonon computations [lOj have treated the spin transi¬ 
tion as in fact a smooth spin crossover, an approach con¬ 
sistent with published experiments. This crossover pro¬ 
ceeds, with increasing pressure, from all Fe ions assuming 
the high-spin state, through to a mixed-spin phase with 
coexisting high-spin and low-spin ions, to eventually all 


ions assuming the low-spin state. 

Previous theoretical work has been based on static 
calculations and has assumed that the mixed-spin state 
is stabilized entropically, yielding a narrow crossover at 
low temperatures that disagrees with experiment. More¬ 
over, experiment and theory have not explored geophys¬ 
ically important properties such as the band structure 
and electrical conductivity. In this Letter, we take a 
different approach that combines first-principles molec¬ 
ular dynamics with free-energy minimization to simu¬ 
late the high-temperature properties of the spin crossover 
directly. Our results reveal a new physical picture of 
the crossover, where the mixed-spin phase is stabilized 
through enthalpy rather than entropy, giving a finite 
broadness for the crossover even at vanishing tempera¬ 
tures. Additionally, we predict the EOS up to the condi¬ 
tions at the base of Earth’s mantle (140 GPa, 4000 K), 
and find that the electrical conductivity of ferropericlase 
reaches semi-metallic values at the bottom of the lower 
mantle, with significant geophysical implications. 

Our simulation setup is built on molecular dynam¬ 
ics (MD) simulations within density functional theory 
(DFT), as implemented in the VASP package [1^. We 
consider a cubic simulation cell of 64 atoms with peri¬ 
odic boundary conditions, adopting an Fe concentration 
of Xpe = 25%, with the Fe ions arranged in a regular su¬ 
perlattice with two nearest-neighbor distances between 
any two neighboring Fe ions. In order to obtain an effi¬ 
cient simulation setup with well-converged values for in¬ 
ternal energy and pressure (within 5 meV/atom and 0.2 
GPa), we sample the Brillouin zone at the Baldereschi 
point [ 2 ^ 211 for a lattice of simple cubic symmetry and 
use a planewave cutoff energy of 500 eV. The projector- 
augmented wave method is used to avoid explicit calcu¬ 
lation of the core electron orbitals. To decide on the best 
feasible approximation to the exchange and correlation 
part of the total energy functional, we compared the EOS 
from conjugate-gradient relaxed static calculations to ex¬ 
periment at 300 K, using the local-density approximation 
(LDA) and two different generalized- grad ient approxima¬ 
tions (GGA), PBE [ 2 ^ and PBEsol [2^. Of these func¬ 
tionals, PBEsol proved clearly superior. 

Unfortunately, the PBEsol functional fails to fully cap¬ 
ture the strong correlation between the 3d electrons of the 
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Fe ions, which is manifested as a spin transition pressure 
of only 18 GPa for xpe = 3.125%, whereas the experi¬ 
mental estimate is closer to 50 GPa [l^. As meta-GGA 
type functionals |24l427l| that we tested brought no al¬ 
leviation to this problem, and as hybrid functional cal¬ 
culations utilizing the exact Fock exchange of the DFT 
system of quasi-electrons are computationally too de¬ 
manding for MD, we use the +U method [1^ to approx¬ 
imate the forementioned correlation effects. Based on 
our calculations on the dependence of the spin transition 
pressure on U, an emp irical estimate for U from optical 
spectroscopy data [l2|, as well as our own hybrid func¬ 
tional [i^ computations, we settled on 17 — J = 2.5 eV. 
We then performed our PBEsol-fl7 MD simulations in 
the NVT ensemble using the Nose-Hoover thermostat. 
Each simulation was run with a timestep of 1.0 fs for a 
total of 10 ps to reach thermal equilibrium, followed by 
10 ps over which all physical time averages were com¬ 
puted. A total of three isotherms, T = 2000, 3000, 
and 4000 K were simulated for compressions that result 
in pressures of approximately 0 to 200 GPa, to encom¬ 
pass existing experimental data and the conditions of the 
lower mantle of the Earth. These dynamic computa¬ 
tions were complemented with a set of static calculations, 
where the crystal structure was relaxed using conjugate 
gradients. 

In order to capture the continuous character of the spin 
crossover and to thus produce a first-principles phase di¬ 
agram of the phenomenon, we minimize the Gibbs free 
energy AG(P, T, /) = G(P, T, /) - G(P, T, 0) at each P 
and T with respect to /, the fraction of Fe ions in the 
high-spin state. As we find a vanishingly small amount 
of intermediate spin (S = 1) Fe in our simulations, 
we define / = where fipe and fipf de¬ 

note the Fe magnetic moment and the same when all 
Fe ions are in the high-spin state, respectively, and () 
denotes an average over Fe ions and time. To map 
AG(P, T,/) = AH — TAS as a function of /, we per¬ 
form constrained-moment and free-moment calculations, 
the former producing a low-spin (/ = 0.0) and high-spin 
(/ = 1.0) phase and the latter producing two mixed- 
spin phases along each isotherm. The enthalpy 77 of a 
given phase / is obtained directly from the MD simu- 
latiqn, as is the electronic contribution to the entropy 
Sei [lo, [sij (we set the electronic temperature equal to 
the ionic temperature). The vibrational entropy Syib and 
entropy Sconf due to site-switching of high spins and 
low spins we evaluate through the method of thermo¬ 
dynamic integration [ 3 ^ . The last contribution to 
the entropy, Smag due to the fully disordered param¬ 
agnetic state of the moments above the Neel tempera¬ 
ture of ^500 K Q, we compute from the expression 351 
Smag = kp + 1)) where Hi is the total magnetic 

moment of Fe ion i, and /cs is the Boltzmann constant. 
We thus obtain AG(P, T, /) at four values of / for each 
P and T, and to find the equilibrium /, we interpolate 
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FIG. 1. Our first-principles phase diagram of the spin 
crossover in ferropericlase. The black line is a geotherm 
from Ref. [3^. 


and minimize AG(P, T, / ) w ith respect to / using a free 
second-order polynomial [3g. 

The resulting phase diagram for the spin crossover is 
presented in Fig. [T] Strikingly, we find a broad pressure 
interval of coexisting high-spin and low-spin ions at all 
temperatures, even at fully static conditions (T = 0 K in 
the phase diagram). Another interesting feature of the 
phase diagram is the weak temperature-dependence of 
the stability field of the mixed-spin phase up to ^3000 K. 
The shape of our phase diagram is thus fundamen tally 
different from previous theoretical work [IS [B [H, 
where the mixed-spin phase was stabilized through an 
ideal mixing entropy, resulting in a completely sharp spin 
transition at T = 0 K. We predict / ~ 0.5 at the core¬ 
mantle boundary, also at odds with previous computa¬ 
tions, which have found significantly smaller high-spin 
fractions. Comparison of our static results to existing 
experimental EOS, XES, and MSB data at 300 K shows 
overall good agreement (Fig. [2^,b). Previous computa¬ 
tions show a much narrower crossover than EOS, XES, 
and MSB data, and our present results. Our results for 
the EOS at all simulated temperatures are presented in 
Fig- Ht: displaying good agreement with experimental 
high-temperature data. 

The finite width of the spin crossover even at vanish¬ 
ing temperatures is due to favorable enthalpy of mixing 
AHmix of the high-spin and low-spin ions (Fig. [3]). We 
trace the favorable AHmix to packing considerations aris¬ 
ing from the volumes of alternating high-spin and low- 
spin Fe-0 octahedra. Due to occupation of Cg orbitals, 
the high-spin octahedron is larger than the low-spin oc¬ 
tahedron, and the Mg-0 octahedron is intermediate in 
size. When high-spin and low-spin Fe in (Mg,Fe)0 are 
brought close together, the system can exploit the will¬ 
ingness of high-spin Fe-0 octahedra to expand and their 
neighboring low-spin octahedra to contract with respect 
to the MgO crystal (see inset of Fig. ED, resulting in lower 
internal energy and forces and hence lower enthalpy than 
expected from ideal mixing of high-spin and low-spin Fe. 
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FIG. 2. a) The EOS of ferropericlase at static conditions compared to experimental data gathered at 300 K [^. [l^ [3^4^ . We 
constrnct the total EOS by interpolating V = V(P, f) linearly between the four spin phases /, where the EOS for each phase 
is a fit to the third order Birch-Murna gha n EOS [42l| . b) Our static result for / compared to XES and MSB data [1, [H, [3 
well as previous computational results [I^ at 300 K. c) Our thermal EOS for all simulated temperatures along with experiment 
at 2000 K [13]. 


These results are consistent with indications that Fe ions 
in (Mg,Fe)0 tend to cluster at high pressure [43|. Our 
finding is, however, in sharp contrast to previous com¬ 
putational work on the spin crossover 2, 18], where the 
high-spin and low-spin ions have been assumed to form an 
ideal solid solution. The favorable AHmix that stabilizes 
the mixed-spin phase in our static simulations persists in 
the dynamical simulations (Fig. [3]). 


Increasing temperature favors the high-spin state be¬ 
cause of the favorable contribution to the free energy 
from the Smag term and to a lesser extent the Syib term. 
Over the range of temperatures that we have considered, 
the magnetic entropy dominates over the electronic en¬ 
tropy, which favors the low-spin state. The mixed-spin 
phase region becomes slightly broader with increasing 
temperature due to the increase in Sconf with increas¬ 
ing temperature. Sconf increases with temperature due 
to the increased occurrence of spontaneous interchanges 
of high-spin and low-spin moments among the Fe sites in 
the mixed-spin phase. 


We find that the vibrational entropy is greater for high- 
spin ions than for low-spin ions. This we relate to the 
shape of the valence charge density of the Fe ion in the 
high-spin state which, considering a sole Fe ion in MgO 
at static conditions and zero pressure, results in a less 
symmetrical Fe-0 octahedron (two axes expanding, one 
contracting) than for the low-spin state (all axes con¬ 
tracting uniformly). This underlying differential distor¬ 
tion, as quantified in our NVT simulations by the dif- 



FIG. 3. Enthalpy (solid lines) and internal energy (dashed 
lines) of mixing of high-spin and low-spin Fe at static condi¬ 
tions, along with enthalpy of mixing at 2000 K and 91 and 
64 GPa (green octagons, left to right, respectively), 3000 K 
and 98 and 73 GPa (orange pentagons), and 4000 K and 106 
and 81 GPa (red squares). Inset: Distribution of octahedral 
volumes at 2000 K and / = 0.50 for Fe-0 octahedra of high- 
spin Fe (blue solid line) and low-spin Fe (red solid line). The 
corresponding static results are shown by the dashed vertical 
lines. The black dash-dotted line denotes octahedral volnme 
for / = 1.0 and 0.0 for static and 2000 K results. 


ference in octahedral quadratic elongation Q between 
the high-spin and low-spin phases, persists at finite tem¬ 
perature, leading to larger mean squared displacements 
of the high-spin Fe ions and hence larger TSmb in the 
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corresponding phase 

The partitioning of Fe between ferropericlase and the 
major lower mantle phase (Mg,Fe)Si 03 perovskite has 
important implications for understanding the structure, 
dynamics, and geochemistry of the Earth’s lower man¬ 
tle 0- We assess the effect of the spin transition on the 
partitioning by computing the ratio \T\{Kf /Ki,q), where 
Kf is the partition coefficient assuming the equilibrium 
/, and Ki,o is the coefficient assuming / = 1.0. Assuming 
no subsequent spin transition in the perovskite, we find 
\n{Kf /Ki o) to lie approximately in the range 0 to 1.5 
along the geotherm [36j, much less than the value of '^10 
estimated by Badro et al. Q. Our much more moderate 
result for the effect of the spin transition on the parti¬ 
tioning appears in better agreement with the relativel 


ely 

(f 


weak pressure-dependence of Kf found in experiment 

The electrical conductivity of the lower mantle is im¬ 
portant for understanding anomalies in Earth’s rotation 
via electromagnetic coupling of mantle and the under¬ 
lying core, and the relationship between observations of 
the geomagnetic field and its source through the filter 
of a potentially conductive mantle. However, no mea¬ 
surements or ab initio predictions of the conductivity of 
ferropericlase at conditions of the deep lower mantle are 
available. Using the Kubo-Greenwood method to com¬ 
pute the electronic component of tr as implemented in 
VASP we find a = 4.0 ± 0.4 x 10^ S/m at con¬ 

ditions close to the bottom of the mantle {P = 136 GPa, 
T = 4000 K), approximately half the recently obtained 
value of 9 X 10“^ S/m for FeO in similar conditions 47j| . 
consistent with the experimental result that a increases 
with Fe concentration [i^. From the electronic density 
of states, it is evident that the metallization of the min¬ 
eral from its initially insulating state is due to the 3d 
electrons of the Fe ions forming broad bands that lead 
to a significant density of states at the Fermi level, an 
effect due to both pressure and temperature. The spin 
crossover itself serves to increase tr, as an increase in the 
concentration of low-spin Fe implies increased density of 
states near the Fermi level of the crystal (3^ . 

The predicted semi-metallic value of electrical conduc¬ 
tivity of ferropericlase at the core-mantle boundary might 
be invoked to explain the highly conductive layer in this 
region inferred from observations of the planet’s nuta¬ 
tions and anomalies therein Assuming the py- 

rolitic volume fraction of 20% for ferropericlase in the 
lower mantle, the presently obtained electrical conductiv¬ 
ity for the mineral, and taking the surrounding perovskite 
phase to be insula ting , the Hashin-Shtrikman minimum- 
maximum bounds ^ for the conductivity of the mixture 
are zero and 5.7 x 10^ S/m, respectively. Adopting the 
maximum value and the half-way value, a simple calcula¬ 
tion shows that respectively 18 or 35 km of lower-mantle 
material is enough to give the required minimum conduc¬ 
tance of 10® S to explain the nutation observations. A 
more highly conductive mantle than previously assumed 


may also require revision of the interpretation of surface 
measurements of the Earth’s magnetic field. 
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SUPPLEMENTAL MATERIAL TO “SPIN CROSSOVER IN FERROPERICLASE 
FROM FIRST-PRINCIPLES MOLECULAR DYNAMICS” 

Details of the free-energy computations 

To compute the full Gibbs free energy for a given isotherm for the high-spin phase (/ = 
1.0), low-spin phase (/ = 0.0), and for / intermediate between 0.0 and 1.0, we start by 
constructing the pressure-volume equation of state (EOS) for each /. For / = 1.0 and 0.0, 
we perform the MD simulations in the canonical NVT ensemble over a range of compressions 
in such a way that the total number of unpaired electrons in the supercell is constrained 
to correspond to either four or zero per Fe ion, respectively. We £t a third-order Birch- 
Murnaghan EOS [l| to the resulting pressure-volume points. 

In order to construct the EOS for / intermediate between 0.0 and 1.0, we hrst perform a 
set of free-moment MD simulations, he., simulations where the number of unpaired electrons 
is not constrained, using various compressions. Each such run produces a certain average 
value for /. We choose two such pressure-volume points, each with a different / intermediate 
between 0.0 and 1.0, to represent a point along the EOS of ferropericlase at each of these 
two /, respectively (see Fig. [1^). We compute the volume for each such intermediate / at 
all other pressures via 

V{PJ) = V{P,0)+g{f)AV{P), (1) 

where AV{P) = V{P, 1) - U(P, 0), and g{f) = [V{P, /) - U(P, 0)]/AU(P) is assumed to be 
constant and equal to the value computed directly from our free-moment simulation for each 
/. The inset of Fig. [1] illustrates that this is a good approximation and a sensible ansatz for 
computing the EOS of mixed-spin phases over the entire pressure range of our calculations. 
Examples of the magnetic structure for each / are shown in Fig. O 

From the EOS for each value of /, we compute the Helmholtz free energy P = E — T{Sei + 
Svib + Sconf) by integrating the thermodynamic identity P = —{dF/dV)T- At this stage, 
the free energy of each spin phase / relative to the other phases is unknown. To resolve 
this, we employ adiabatic switching within the Kirkwood coupling scheme to compute the 
free energy of a given phase / relative to / = 0.0. The joint potential energy function for 
switching from / = 0.0 to an arbitrary value of / with the coupling parameter A going from 
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FIG. 1. The final result for the EOS at each value of /, this example being for T = 3000 K. 
The EOS for / = 1.0 and / = 0.0 are htted to data points from simulations with constrained 
moments (marked by circles). The EOS for / = 0.65 and / = 0.47 are constructed starting from 
free-moment simulation points (marked by crosses). Inset: Explicit test of the method used to 
construct the EOS for / intermediate between 0.0 and 1.0, at static conditions. The large symbols 
indicate points from which the EOS were constructed, small symbols indicate points calculated 
explicitly after the construction of the EOS. 


0 to 1 is 

E{X) = Ef=o,o + X{Ef — Ef=o,o). 


( 2 ) 


The free energy difference between / = 0.0 and a given / is computed as the integral 


AE — Ef — F/=o.o — [ dX{Ef — £'/=o.o)a, 

JO 


(3) 


For / = 1.0, we use A = 0 and 1 to constrain the integrand at the end points. Additionally, 
we employ the fluctuation term derived from thermodynamic perturbation theory 

((F'/=i.o - F'/=o.o)^)a=i - {Ef=i.o - Ef=o.o)l=i 




(4) 


which, from simple geometrical arguments, is equal to —1/2 x {d{Ef=i,o — £'/=o.o)/<^A)a=i, 
provided the perturbation term of Eq. IHis small (we hnd 4 to 13 meV/atom throughout all 
calculations). We thus constrain the slope of the integrand at A = 1, and using the three 
constraints, take a fully constrained second-order polynomial to model the integrand. For / 
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FIG. 2. Magnetic moment of each Fe ion in the 64-ion supercell as a function of time after the 
initial 10 ps equilibration, in units of the high-spin Fe moment. Data are shown for / = 0.0 (a), 
/ = 1.0 (b), / = 0.76 (c), and / = 0.49 (d) at T = 4000 K. The inset in (d) is the corresponding 
data for / = 0.39 at T = 2000 K. 


intermediate between 0.0 and 1.0, we use A = 1 and the fluctuation term of Eq. 01 as in 
In Fig. [3l we present an example of the integrand and an assessment of the relative precision 
of the approximations detailed above, showing that overall the method looks sufficiently 
precise for the present purposes. 

Finally, we add on the magnetic entropy Smag to the free energy of each spin phase /, 
and move on to analyze the results in the NPT ensemble, i.e., in terms of Gibbs free energy. 
To hnd the equilibrium high-spin fraction at each pressure and temperature, we interpolate 
AG(F, T, /) = G{P, T, f) — G{P, T, 0) with a free second-order polynomial between the four 
values of / and hnd the minimum of that polynomial in the interval / G [0,1] (See Fig. 0] 
for an example). We consider this method of interpolation a more neutral choice than using 
a third-order polynomial or natural splines, both of which seemed to produce spurious and 
artihcial features into the interpolation. 
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FIG. 3. a) The integrand of Eq. [3] at T = 3000 K at a lattice constant of 3.93 A for / = 1.0, 
using the approximation of constraining the integrand only at A = 0 and 1 (linear fit) and the 
approximation of constraining additionally the slope of the integrand at A = 1 (quadratic fit). The 
values of the integrand at A = 0 and 1 are shown for reference, b) The value of the integral of Eq. [3] 
for / = 1.0 for all isotherms and different levels of approximation. Using only A = 0 or 1 results 
in overestimation or underestimation of AT, respectively. The linear and quadratic fit agree well 
with the approximation of using A = 1 and the fluctuation term of Eq. m implying good precision 
for the chosen approximations. 
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FIG. 4. Finding the equilibrium value of / at T = 3000 K and P = 120 GPa through second-order 
polynomial interpolation and minimization. 


Explaining the vibrational entropy results 


Populating the spin-up and spin-down t 2 g and Cg states of a single Fe ion within an MgO 
lattice unevenly to produce the high-spin conhguration entails a valence electron density 
around the ion which results in a non-uniformly distorted Fe-0 octahedron, contrary to the 
situation in the low-spin state. We quantify the distortion of the Fe-0 octahedra in the 
dynamic simulations by using the mean octahedral quadratic elongation (Xoct) jd], and hnd 
Fe-0 octahedra at / = 1.0 to be generally more distorted than their counterparts at / = 0.0. 
In Fig. [5] we show that the difference in mean squared displacements of the Fe ions between 
/ = 1.0 and 0.0 are correlated with the difference in (Xoct) between the two phases. This 
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FIG. 5. Explaining the results for the vibrational entropy difference between / = 1.0 and / = 0.0, 
this example being at T = 3000 K. The difference between the two phases in mean squared 
displacements (r(t)^) of the Fe ions with respect to their mean positions is correlated with the 
difference in octahedral elongation, ie., the degree of distortion of the Fe-0 octahedra. The 
difference in vibrational entropy TS^n, between / = 1.0 and / = 0.0 is of the order of 15 meV/atom 
here. The mean octahedral volume is the same for both phases in these NVT simulations. 

implies that the ionic network at / = 1.0, he., with all Fe ions in the high-spin state, is more 
distorted than at / = 0.0, with all Fe ions in the low-spin state, leading to larger nuclear 
excursions and hence the observed larger vibrational entropy. 


Results on Fe partitioning 

Results for the ratio of the two partition coefficients described in the main article are 
shown in Figs|6^ and b. 


Electronic density of states and metallization 

In Fig. [TJ we show the electronic density of states (DOS) of ferropericlase. Increasing 
temperature and pressure causes the DOS near the Fermi level to increase significantly, and 
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FIG. 6. a) Ratio of partition coefficient of Fe between ferropericlase and perovskite assuming the 
equilibrium high-spin fraction (Kf) and all Fe in the high-spin state (iFi.o). b) Bilinear interpolation 
of the results of a), along with a geotherm ^ (white solid line). 


hence electrical conductivity to increase, eventually up to semi-metallic values. Collapsing 
the Fe moments to the low-spin state, as during the spin crossover, further increases the 
DOS near the Fermi level. 
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FIG. 7. Total and projected DOS of ferropericlase illustrating the effect of pressure and temperature 
as well as high-spin fraction on the DOS near the Fermi level. The plots are shifted so that the 
Fermi level is at 0 eV. 
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